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We have studied a model of a random walk in a quenched random environment. In addition to 
featuring anomalous diffusion and localization, for special regimes of disorder parameters the particle 
density decomposes into multi-Gaussian structure while its cumulative distribution is normal. We 
explain the observed fine structure of the density and point out its significance to experiments. 
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I. INTRODUCTION 

Diffusion is ubiquitous in nature. The first theoret- 
ical description of diffusion processess such as Brown- 
ian motion of particles in colloidal suspensions was de- 
veloped by Einstein and Smoluchovski [TJ [2]. Normal 
particle diffusion is characterized by the fact that the 
mean square displacement of the corpuscle grows linearly 
with time and the cumulative probability distribution is 
asymptotically that of a Gaussian density as dictated 
by the Central Limit Theorem (CLT). However, under 
special circumstances the situation may change dramati- 
cally, e.g., in Sinai's model the growth of the variance be- 
comes logarithmically slow [3] , and the particles may be- 
have super- or subdiffusively or even become localized — a 
phenomenon predicted by Anderson for electron waves in 
disordered crystals [4 and most recently observed with 
Bose-Einstein condensed matter waves in a controlled 
disorder [5j [6] . 

Over the past few decades, the study of wave and par- 
ticle propagation in random media has found a broad 
range of applications in diverse fields including medicine, 
optics, materials research, quantitative finance, and the 
biology of epidemics and genetic selection [TJ [8] . The ap- 
proach to understanding random media has often been 
through random walks in random environments (RWRE) 
[3] [9tiTTj . In these models the environment describes the 
local propagation laws, which model the local properties 
of the inhomogeneous medium in question [T2HT4] . In a 
quenched (frozen) environment, the degree of stochastic 
freedom in the study of an RWRE is limited rendering 
the mathematical analysis of the problem challenging. 
In this setting, major advances have been obtained only 
very recently [T5HT7] . 

In this paper, we consider combined ballistic and dif- 
fusive particle dynamics [18] in a spatially random but 
temporally fixed environment. Our model originates in 
the study of certain extended dynamical (i.e., deter- 
ministic) systems obtained by coupling together an in- 
finite sequence of simple chaotic dynamical systems [19] . 
With random initial conditions these can be reduced to 
RWREs in quenched environments by introducing a so- 
called Markov partition of the phase space. The model 
in question is also related to quenched compositions of 
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FIG. 1. (Color online) Schematic of our model of one- 
dimensional random walk in a quenched random environment. 
The labeled rectangles denote the sites of the discrete array 
and the arrows illustrate the site to site jump probabilities. 
Sites with the same blue or green color (dark or light gray) 
have the same fixed transition probabilities. The notation is 
explained in the text. 



chaotic maps [20] . 

A natural question regarding such a stochastic process 
is whether the resulting particle distribution is asymptot- 
ically Gaussian. By virtue of Sinai's example [3 , 21, 22 
one may not always expect an answer in the affirma- 
tive, even when the environment is not quenched. In our 
model the parameter space splits into subregions corre- 
sponding to different qualitative behavior; in some re- 
gions with ballistic behavior or parametric symmetries 
the cumulative distribution is normal, while in others it 
is not. Even in the former case the particle density is 
not normal, and under certain circumstances features a 
multi-Gaussian structure — a striking phenomenon first 
observed in [19]. Here we also explain the source of this 
observation. 



II. MODEL 
A. Setup 

We consider one-dimensional RWRE as illustrated in 
Figure [1] Within a single time step, from a given site 
there are prescribed probabilities of jumping into either 
of the two nearest neighbor sites. It is also permissi- 
ble to remain in the same site with nonzero probability. 
Put more precisely, let X n E Z = {. . . , —1, 0, 1, ... } be 
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the position of the particle at time n G N = {0, 1, ... } 
and assume Xq = 0. The particle jumps about the 
lattice Z according to Markov transition probabilities 
Pij = P(X n+ i = j\X n — i). To describe these tran- 
sition probabilities completely, we assign to each site i a 
triplet uoi = (qi,ri,pi), where n and pi are the proba- 
bilities of jumping to the left, staying put, and jumping 
to the right, respectively, if the current site is i. That 
is, q { = Pij-i, r» = p iti , Pi = p*,t+i. In our model 
Qi + r i + Pi = 1? because the length of a jump cannot 
exceed one. We will also assume that qi > and Pi > 0. 
The collection uj = (. . . , gj_i, cjo, cji, . . .) is called an en- 
vironment. 

Let {(q £ ,r £ ,p £ ) : £ G A} be the collection of all pos- 
sible triplets that each uji may be equal to. Then, the 
triplets uji are drawn as i.i.d. trials from this collection 
with equal probability. Hence, for each site z, we have 
uoi = {q l ,r l ,p l ) for some label £ G A, and this choice is 
independent of all the other sites. Once every uji has been 
randomly picked, it is frozen for good. An environment 
uj — (. . . , cj_i, c^o, cji, . . .) generated in this way is called 
a quenched random environment. 

For example, we can take A = {A, B} if there are only 
two different kinds of triplets in the environment. Having 
chosen uji = (q £ ,r £ ,p £ ) for a site z, where £ is either A or 
£?, we call interchangeably uji and £ the label of the site. 
We also say that the site z is of type £. 

The sequence (X n ) n > of random variables is a Markov 
Chain and describes a random walk in a quenched ran- 
dom environment. We are interested in the probabil- 
ity distribution of Z n = (X n — EX n ) / \/ Var X n for typ- 
ical environments uj. In particular, CLT is said to hold, 
if lim n ^ 00 P(Z n < t) = $(£) for all t G R. Here 
<$>(t) = -A= f* e _s / 2 is the cumulative distribu- 
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tion function of the standard normal distribution A/"(0, 1). 
However, our focus is on the form of the density function 
P(X n = k) itself, which admits surprising features. 

B. Implementation 

To compute the probability distribution of X n we have 
employed two complementary numerical methods. In the 
first method we calculate the exact probability distribu- 
tion of X n by computing the nth power of the transi- 
tion probability matrix of the Markov Chain. Since the 
lengths of the jumps are at most one and the walks start 
from Xo = 0, we have \X n \ < n. Therefore, it suffices 
to compute powers of a finite matrix to find out the dis- 
tribution of X n for any fixed n. Depending on the cho- 
sen parameters the distribution may in fact spread quite 
slowly with n and a relatively small matrix can be used 
if n is not too large. To study larger values of n, we use 
a Monte Carlo algorithm. Starting at and proceeding 
recursively, following a position i = X n the next position 
X n +i G{z — l,z,z + 1} is picked randomly according to 
the probabilities The procedure is repeated 

and averaged over many independent particle trajecto- 



ries to obtain an approximation for the distribution of 
X n . In the figures presented the former exact method 
has been applied. 



III. PRELIMINARIES 

A. Recurrence, transience, and ballisticity 

The walk is called recurrent if, with probability 1 (w.p. 
1), it returns to 0. In other words, at some later (ran- 
dom) time n > 0, X n = 0. This clearly implies that, 
in fact, the walk returns to infinitely often. For our 
model it follows from [10 that recurrence is equivalent 
to lim sup^^^ X n = oo and lim inf V^oo X n = — oo (w.p. 
1). The walk is called transient if it is not recurrent. In 
this case one of only two things may happen [10] . Either 
lim n ^oo X n -oo (w.p. 1) or lim^^ X n = oo (w.p. 1). 
In the former case the walk is transient to the left, in the 
latter it is transient to the right. While a transient walk 
does escape to infinity, it may do so extremely slowly, so 
that ^ — >• 0. A transient walk is called ballistic if the 

n 

limit linin^oo ^ exists and is nonzero. 

The quantity // = % is called the bias of the label 
£ G A. The recurrence and transience properties of the 
walk can be conveniently described with the aid of the 
biases of the labels present in the environment; see below. 



B. Reversible state tt 

There is a natural potential function V associated 
with the environment uj. Indeed, SVi = log can be 
thought of as a potential difference between the sites i 
and i + 1: if SVi > 0, the likelihood of jumping from 
i to i + 1 is less than vice versa, and the potential at 
i + 1 should be higher than at i. To fix an arbitrary con- 
stant, we set Vo = 0. We then define V^+i = V{ + SVi 
for all z, such that Vi = X^o<j<i ^3 f° r ^ > and 
Vi = — Xli<j<o f° r ^ < 0- A similar potential was 
considered in [3]. 

The Gibbs-like state tt defined by 7r^ = e~ Vi satisfies 
the detailed balance condition TTiPij = 7tjPj,i (|z — j\ = 
1) and is therefore a reversible stationary state for the 
Markov Chain. This follows immediately from 

e -V z+1 = e -( Vi +6Vi) = e -Vi gM±l _ 

Notice that tt cannot be normalized to a probability dis- 
tribution, because in general Y^i^i = 00 ■ I n ^ ac t? the 
potential Vi may assume arbitrarily large negative val- 
ues sufficiently far away from i = 0, which makes iXi 
unbounded. As we shall see, we can nevertheless learn 
about the distribution of the walk from the state tt. It 
is particularly useful in situations in which the walk is 
recurrent. 
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C. Stationary state li 

Now suppose the walk is transient to the right. Be- 
cause of this property, it will visit every site i finitely 
many times. Let Ni be the total amount of time spent at 
site i > 0. Notice that, for every realization of the walk, 
there is precisely one more jump from site i to i + 1 than 
vice versa. This follows from transience to the right and 
the fact that the maximum jump length is one: every 
jump from i + 1 to i must sooner or later be followed by 
a jump from i to i + 1. On the other hand, the expected 
number of jumps from i to i + 1 is E(Ni)pi and the ex- 
pected number of jumps from i + 1 to i is E(i\^ + i)g^ + i. 
Thus, we have the relation E(A^)^ — 1 + E(7V^i)<^+i. 

The above reasoning motivates the following con- 
struction. Demand that \±i = ^-(1 + /i^+i^+i) for 
every i G Z. Writing p i = ^ for brevity, the expression 

to = ^( 1 + Pi+l+Pi+lPi+2+Pi+lPi+2pi+3 + . • •) follows. 

Transience and ergodicity of the i.i.d. environment guar- 
antee that the latter series converges. It is easily checked 
that the state \i so defined is a stationary state: the 
global balance condition pi-iTTi-i + r^i + g^ + i7r^ + i = 7^ 
holds for every i. 



IV. RESULTS 
A. Anomalous probability density 

In the case where the bias // has the same value for 
all labels £ G A, we find the probability distribution to 
decompose into |A| distinct Gaussian functions, where 
| A| is the number of elements in the set A. The equal bias 
condition can be satisfied only if r £ is different for every £, 
which means that staying put must be allowed if |A| > 1. 
Figure [2] shows the scaled and centered probability 
densities of X n for the equal bias cases in which the sets 
of possible labels uji are (a) {( ^, ^, 
( D ) {( A' iu' A)' ( A' in' ( A' A' M)} and (c) 

Jf^_ 16 _2n /_3_ 14 _3_n 12 __4_ \ , _6_ _8_ _6_ n / 10 _0_ 10 J 
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In each frame, a different marker is assigned according 
to the environment label ujx n at the endpoint X n . 
Joining the points of same marker type would result in 
perfect Gaussian curves which collapse on top of each 
other when multiplied by a constant. The ordinates 
of the horizontal lines, drawn to guide the eye, follow 
the same ratio as the right transition probabilities p l 
(equivalent ly the left transition probabilities q £ ) of the 
labels on the corresponding curves. The latter curves 
themselves are similar; in each frame, multiplying any 
such curve with a certain constant results in perfect 
overlap with the solid Gaussian shown. Notice that 
in the equal bias case pi = const • Thus, at the 

Pi 

local level, the particle distribution closely reflects the 
structure of the stationary state \i. 

This picture is stable under variation of n: there are 
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FIG. 2. (Color online) Probability densities of Z n , where 
n — 10 5 , for a symmetric (a) two (blue and black dots) (b) 
three (blue, black and red dots) and (c) five (blue, black, red, 
green and cyan dots) label system, where the color (shade of 
gray) specifies the type of the endpoint label ux n ■ The tran- 
sition probabilites are listed in the text. The ordinates of the 
horizontal lines, drawn to guide the eye, obey the same ratios 
as the right transition probabilities pi of the corresponding 
site types £ G A. The solid (magenta) curve in each frame is 
a Gaussian function with zero mean and a variance measured 
from the data. Each dotted curve coincides with the exact 
Gaussian when multiplied by a proper constant; see text. 

rapid oscillations from site to site, so the centered, scaled, 
and interpolated probability density does not converge 
to any function. Instead it has several smooth Gaussian 
"components" . Altering the relative concentration of the 
site types only changes the concentration of points on 
each Gaussian. We emphasize that the multi-Gaussian 
structure only appears in the presence of the fixed envi- 
ronment. If one averages the distributions over many en- 
vironments (in addition to all possible trajectories which 
are already included), the data collapses to a single Gaus- 
sian function. The multiple Gaussian curves persist as 
long as the biases // are equal. From here on we focus on 
the two label case A = {A, B}. Increasing the number of 
site types does not affect our main conclusions. 

B. Analysis of observed densities 

The behavior of the probability density may be under- 
stood intuitively by considering the potential function 
Vi. Locally, particles tend to accumulate in the valleys 
of the potential function and correspondingly potential 
peaks repel particles. The bias pi yields the relative ten- 
dency of the site i to transport particles to the left or 
to the right once they leave the site. Domain bound- 
aries where neighboring sites tend to transport particles 
in the opposite directions correspond to the local extrema 
of the potential function. Knowledge of the structure of 
the potential, which is straightforward to obtain once 
the environment is specified, turns out to be enough to 
infer qualitative features of the density such as multi- 
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position 



FIG. 3. (Color online) Probability densities for (a) p A — p B ^ 
1, (b) 1 > p A > p B and (c) p A p B = 1 / p A . In (a) potential 
function, the stationary state p, and the density are plotted 
using different colors (shades of gray) and are tagged next to 
each of the curves. In (b) the density and a Gaussian curve 
modulated by the stationary state p closely follow each other 
(potential is not visible in this frame). Frame (c) shows the 
potential, density, and the reversible state tt as indicated next 
to each of the curves. Horizontal lines are drawn to guide the 
eye. 



Gaussianity without having to calculate or measure the 
density itself. 

Figure [3^a) shows a probability density together with 
the potential function for an equal bias case, p A — p B ^ 
1. In addition, the stationary state p is plotted. The po- 
tential forms two linear potential levels resulting in the 
observed double Gaussian density. For general values of 
the biases, certain configurations of the environment la- 
bels yield clearly distinguishable linear potential levels. 
Correspondingly, the particle density has equally many 
Gaussian envelopes. For example, if the environment is 
made periodic, the potential (and hence also the den- 
sity) always has multi- level structure also when p A ^ p B . 
This suggests to us that, whenever the CLT holds, the 
density in fact is a modulated Gaussian. The modulat- 
ing factor at each site is a function of the environment 
but not of time. Although the probability of sampling 
a perfectly periodic environment is 0, such environments 
may be physically important since many structures in 
nature exhibit approximate periodicity, including crystal 
structures in various solid state systems. Periodic envi- 
ronments can also be engineered on purpose. 

The evidence points to the stationary state pi being 
the modulator in the ballistic CLT regime. Indeed, the 
density in Figure [3^a) is perfectly reproduced (less nor- 
malization) by multiplying p with a Gaussian having the 
measured mean and variance. Further, Figure [3^b) shows 
such modulated Gaussian which closely follows the actual 
density for a more generic transient case. The potential 
function is not visible in Fig. |3|b) since it falls outside 
the scale shown in the figure. 




FIG. 4. (Color online) Schematic diagram displaying param- 
eter regions exhibiting qualitatively different system behav- 
ior (assuming r £ = 0). On the blue (light gray) hyperbola 
p A p B — 1 the walk is recurrent and shows Sinai diffusion, 
except precisely at p A — p B = 1 where CLT holds. In re- 
gion (I) the walk is transient to the left. In regions (II) and 
(III) the walk is transient to the right. Inside the black (dark 
gray) circle (pa) 2 + (pb) 2 = 2 the CLT holds. In the region 
(III) between the hyperbola and the circle, the CLT does not 
hold. When r £ ^ 0, double- Gaussian densities emerge on the 
magenta line p A — p B . 

In Figure [S^c), p A p B = 1 / / and the particles be- 
come localized in the potential wells and make occasional 
quick journeys from one well to the next, which is anal- 
ogous to tunneling through potential barriers. In this 
recurrent case, an excellent prediction for the density is 
given by the reversible state ir as illustrated in Fig[3^c). 
Notice, in particular, how the valleys and peaks in the 
potential are correlated with the maxima and minima, 
respectively, of the probability density. 



C. Phase diagram 

Let us collate a few known results relevant to our two 
label model, A = {A, B}. With the exception of re- 
currence, transience, and ballisticity, these results have 
been proved for the case in which staying still is forbid- 
den (r £ = 0). Recall that p £ = ^ is the bias of the site 
type l = A,B. 

According to [10j[T4], if p A p B — 1, the walk is recur- 
rent. If p A p B < 1, the walk is transient to the right. If in 
addition p A -\-p B < 2, the walk is ballistic. By symmetry, 
left transience corresponds to p A p B > 1. Below we will 
discuss recurrent and right transient cases. 

In the recurrent case, it has been shown [23] that the 
CLT holds under the symmetry p l = 1 (£ = A, B). The 
reason is that the stochastic process (X n ) is a martingale 
due to the symmetry. In the absence of this symmetry the 
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variance grows according to the Sinai law (logn) 4 . Sinai 
[3] studied walks for which the probabilities p i = \ — q i 
(ri = 0) are i.i.d. random variables whose values are not 
quenched, meaning that the walks are averaged over all 
environments. His result is that, in the recurrent case the 
variance of X n grows as (logn) 4 , which is astoundingly 
slow compared with the usual n. The limit distribution is 
not Gaussian: it has been computed explicitly by Kesten 
and Golosov |2TJ[22]. However, if one allows > and 
imposes the additional inter-site symmetry pi = q i+ i for 
all z, which yields a flat potential (SVi — 0), the usual 
CLT does hold [24]. 

Under the assumption r £ = and p A p B < 1 (tran- 
sience to the right), Goldsheid [17 has shown that 
(p A ) 2 + (p B ) 2 < 2 suffices for the CLT. Now, if p A > 1, 
then p B < 1 by transience, such that the equation 
(p A ) s + (p B ) s — 2 has precisely one positive solution 
s > 0. If < s < 2, the statistical behavior of X n is com- 
plicated and there is no limit distribution [25]. Indepen- 
dently of Goldsheid, Peterson [25] has also obtained the 
CLT, under stronger assumptions which include s > 2. 

Figure [4] summarizes the above knowledge of the case 
r £ = {I = A, B) in the form of a schematic phase di- 
agram of the system in terms of the biases p A and p B . 
Superimposed is the double- Gaussian region we have ob- 
served when r £ > 0. 

For p A = p B , double- Gaussian density functions 
emerge. If p A w p B , both the probability density and 
the underlying potential exhibit a double structure. How- 
ever, only for p A = p B are the two densities strictly Gaus- 
sian. The larger the difference in the biases the larger 
the local fluctuations in the components of the probabil- 
ity density such that eventually the fluctuations grow to 
become of the order of the amplitude of the curves them- 
selves and consequently the double Gaussian structure is 
washed away as illustrated in Figure |5|a), where uji G 

{(li h ra). When p a < 1 and p b ^ !' or 

vice versa, irregularity of the distribution even at the cu- 
mulative level is observed in simulations, as shown in Fig- 
urep|b) , where uji G { ( ^ , ^ , ^) , ( ^ , ^ , ^ )}. Whether 
this irregularity eventually disappears or persists for- 
ever depends on whether the parameters are in the CLT 
regime — which at least for = is (p A ) 2 + (p B ) 2 < 2 
— or not. 



D. Scaling of variance 



Figure [6] shows the variance of X n as a func- 
tion of n computed for three different parameter 

sets- m u- G U^- -^-) ^ -^-)\ (ii) 

bt^ub. ^1; w % t \ V 100 ' 100 ' 100 / ' v 900 ' 900 ' 900 / / ' W 

ri f= {( 1 89 10 ^ ( 1 889 -10- U anrl (\\\\ , > d 
w% ^ I v 100 ' 100 ' 100 / ' v 900 ' 900 ' 900 / J ' clliU v m / w% ^ 

{(iM> OT' Iuo)' (9M' §§§' 9^)}- F° r W tne Potential 
function has horizontal two-level structure and the vari- 
ance grows linearly with n exhibiting normal diffusion. 

In case (ii) the two-level potential is tilted induc- 
ing global drift to the system and the resulting vari- 
ance shows prominent persistent oscillations. The sys- 




FIG. 5. (Color online) Cumulative distribution function for 
(a) 1 > p A > p B and (b) p A p B = 1 / p A . Insets show 
the corresponding probability densities where the solid curves 
are Gaussian functions with zero mean and whose variances 
correspond to the data. For both frames n = 10 5 . 



tern mimics normal diffusion on average but alternates 
between sub and superdiffusive behavior. This poten- 
tially transient behavior is understood to be caused by 
the probability mass propagating through focusing and 
defocusing regions of the fixed potential landscape. For 
practical purposes the fluctuations are quite prominent 
and should be significant for experiments. 

The third case (hi) demonstrates the effect of strongly 
trapping quenched environment whence the particles 
need to tunnel through arbitrarily high potential barriers 
and the system exhibits localization, as expected on the 
basis of the literature on the case r £ = 0. By the same 
token, we typically observe localization of the density of 
the kind displayed in Figure pi when p A p B w 1 (excluding 

fVio m Q-rfin (TqIo pqco — /i — "|\ . 



the martingale case p A = p B ~= 1), although by abstract 
results on the r £ = model this is likely to be a transient 
phenomenon if (p A ) 2 + (p B ) 2 < 2 (region II in the phase 
diagram) . 



V. DISCUSSION 

We have investigated a model of a random walk in a 
quenched random environment. In the existing litera- 
ture, particle distributions have been studied at a coarse 
diffusive scale. Indeed, the CLT gives asymptotic approx- 
imations of the probabilities P(a v /n < X n — EX n < by/n) 
for fixed pairs a < b in the large n limit. In our work we 
have investigated substantially finer details of the distri- 
bution by studying P(X n = fc), for fixed values of k. 

It turns out that the latter local distribution has an un- 
expectedly rich, highly oscillatory, nature. This is true 
even in the CLT regime, in spite of the fact that the 
cumulative distribution obtained by integrating the par- 
ticle density over large length scales reduces to a perfectly 
Gaussian cumulative function. Such prominent and per- 
sistent oscillations are significant to experiments in which 



6 



10 4 













<^(logn) 4 




, ,/<\ 





10- 1 — — — — ^ — — — ^ — — I 

10° 10 1 10 2 10 3 10 4 10 5 



FIG. 6. (Color online) Variances as functions of time for (i) 
p A = p B = 1, (ii) p A = p B / 1, and (iii) p A p s = 1 / p A . The 
functional forms of the solid curves are marked in the figure. 
With nearly equal bias the variance grows linearly with time 
exhibiting normal diffusion. When at least one of the biases 
is different from unity the variance begins to fluctuate. In 
the strongly trapping case (opposite bias) particles become 
strongly subdiffusive and localize. 



information regarding the scaling limit of the system is 
not sufficient. In this paper we have obtained theoretical 
understanding of them. 

For equal-bias situations, as well as for periodic en- 
vironments with generic biases, the particle density is 



observed to break up into multiple Gaussian envelopes 
whose relative amplitudes are given by the ratio p £ /p £ 
for labels G A. This ratio is found to be correctly 
predicted by the stationary state /i. Generically, for tran- 
sient cases the density can be predicted by multiplying \i 
with the Gaussian whose mean and variance match those 
of the observed distribution. In recurrent cases we find, 
as one would expect based on related literature, that the 
reversible state tt provides the most accurate prediction 
for the particle density. 

For the sake of concreteness, we have used traditional 
terminology in that the walk represents a particle in a 
medium. This should not be construed as a restriction 
of applicability to other types of processes. Nor should 
it be mistaken to mean that the basic constituents, the 
sites (or 'cells'), must be small let alone microscopic. We 
believe that the observations made here can be extended 
to a range of processes occurring in random or otherwise 
disordered environments. In this paper, we have stud- 
ied an ensemble of noninter acting particles. In future, 
it will be interesting to investigate the universality, per- 
sistence and stability under weak particle interactions of 
the observed multi-Gaussian densities in both classical 
and quantum systems in the presence of quenched envi- 
ronments. 
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